mypal = c("#1f77b4", # ast
          "#ff7f0e", # end
          "#2ca02c", # ext
          "#d62728", # inh
          "#9467bd", # mic
          "#8c564b", # opc
          "#e377c2"  # oli
)
mypal = setNames(mypal,c("ast","end","ext","inh","mic","opc","oli"))

cell_names = c("Astrocytes", "Endothelial cells", "Excitatory Neurons", "Inhibitory Neurons", "Microglia", "Oligodendrocytes", "OPCs")
cell_names = setNames(cell_names, c("ast","end","ext","inh","mic","oli","opc"))

pheno_list = c("cogng_demog_slope"="gaussian", # Cognitive decline slope. Remove the effect of demog 
               "cogng_path_slope"="gaussian", # Resilience, removed the effect of path + demog
               "tangles_sqrt"="gaussian", # Tangle density - Mean of 8 brain regions
               "amyloid_sqrt"="gaussian", # Overall amyloid level - Mean of 8 brain regions
               "gpath"="gaussian", # Global burden of AD pathology based on 5 regions
               "tdp_cs_6reg"="gaussian", # TDP-43, 6 region severity summary
               "ad_dementia_status"="binomial" # Clinical AD # CT = MCI + NCI 
               )

Load sn

sn_columbia <- read.table(paste0(ids_dir, "projid_sn_columbia.txt"), header = T)
sn_columbia$projid <- sprintf("%08d", sn_columbia$projid)

load(paste0(work_dir,"/MIT_pseudobulk_norm.RData"))

Read sn modules

se_results_dir = "/pastel/projects/speakeasy_dlpfc/SpeakEasy_singlenuclei/2nd_pass/snakemake-sn/results/"
# Load sn modules
sn_modules = data.frame()
for(cell_type in c("ext","inh","ast","oli","mic","end","opc")){
  # Filter modules with lest than 30 genes
  modules_file = read.table(paste0(se_results_dir,cell_type,"/geneBycluster.txt"), header = T)
  modules_file = modules_file[, c("ensembl", "gene_name", "cluster_lv3")]
  modules_file$module_clusters = paste0(cell_type,"_M", modules_file$cluster_lv3)
  modules_file$cluster_lv3 = NULL
  
  modules_file = modules_file %>% dplyr::group_by(module_clusters) %>%
    dplyr::mutate(module_size = n()) %>% dplyr::filter(module_size >= 30) %>%
    dplyr::mutate(cell_type = cell_type)
  
  sn_modules = rbind(sn_modules, modules_file)
}
# length(unique(sn_modules$module_clusters)) # 193
# modules by cell type
table(unique(sn_modules[,c("module_clusters","cell_type")])$cell_type) # 7
## 
## ast end ext inh mic oli opc 
##  26  26  29  24  30  30  28

Project sn modules in the MIT data

We already have the modules from the Columbia sn-ROSMAP, so I’ll use them as a gene list and calculate the eigengenes and module average expression for the MIT sn-ROSMAP dataset. Then, we’ll run regressions with AD-related covariates to check for associations.

Heatmaps: Bonferroni adjusted by column.

Astrocytes

We keep only the MIT samples NOT in Columbia sn.

# Astrocytes 
net_output = unique(sn_modules[sn_modules$cell_type == "ast", ])
expr_mtx_sel = as.data.frame(celltype_exp$ast$tmm_voom)[, ! colnames(celltype_exp$ast$tmm_voom) %in% sn_columbia$projid]
rownames(expr_mtx_sel) <- gsub("(*.)\\.(.*)","\\1", rownames(expr_mtx_sel))
expr_mtx_sel$ensembl = rownames(expr_mtx_sel)
gene_mod_map = data.frame(ensembl = rownames(expr_mtx_sel))
gene_mod_map = gene_mod_map %>% inner_join(unique(net_output[,c("ensembl","module_clusters")]), by = c("ensembl"))
gene_mod_map$cluster_lv3 = gsub("(.*)_M(.*)","\\2",gene_mod_map$module_clusters)

expr_mtx_sel_ord = expr_mtx_sel[match(gene_mod_map$ensembl,rownames(expr_mtx_sel)),] # order the matrix
identical(rownames(expr_mtx_sel_ord), gene_mod_map$ensembl) # must be TRUE 
## [1] TRUE
colors_mod = gene_mod_map$cluster_lv3
expr_mtx_sel_ord$ensembl = NULL
expr_mtx_sel_ord_t <- data.frame(t(as.matrix(expr_mtx_sel_ord)))
external_ME = moduleEigengenes(expr_mtx_sel_ord_t, colors = colors_mod)
mod_average = external_ME$averageExpr # data.frame

# Save results 
save(external_ME, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleEigengenes_lv3_ast.Rdata"))
write.table(external_ME$eigengenes, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleEigengenes_lv3_ast.txt"), sep = "\t", quote = F, row.names = T)
write.table(external_ME$averageExpr, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleaverageExpr_lv3_ast.txt"), sep = "\t", quote = F, row.names = T)

Regressions

data4linear_reg <- mod_average # match codes name 
phenotype_dt = phenotypes[match(rownames(data4linear_reg), phenotypes$projid), ]

all(rownames(data4linear_reg) == phenotype_dt$projid) # Must be TRUE. Match IDs
## [1] TRUE
# hist(phenotype_dt$amyloid_sqrt)
message(paste0("Unique participants: ", length(unique(phenotype_dt$projid)))) 
## Unique participants: 191
res_test = run_module_trait_association(data4linear_reg = data4linear_reg, 
                                        phenotype_dt = phenotype_dt, 
                                        pheno_list = pheno_list, 
                                        covariates = c("age_death","msex", "educ")) # covariates to adjust
matrix_rsquared = res_test$matrix_rsquared
matrix_pvalue = res_test$matrix_pvalue
save(res_test, file = paste0(work_dir, "replication_analysis_filt/MIT_res_test_ast.Rdata"))

plot_module_trait_association_heatmap(res_test)

Top results

Top result by covariate.

# res_test$all_stats_df 
createDT(res_test$all_stats_df %>% group_by(phenotype) %>% slice_head(n = 1))

Nominal pvalue

createDT(matrix_pvalue)

Significant results

Threshold: At least one module with adjusted pvalue < 0.05.

plot_module_trait_association_heatmap(res_test, show_only_significant = T, signif_cutoff = 0.05)

Microglia

net_output = unique(sn_modules[sn_modules$cell_type == "mic", ])
expr_mtx_sel = celltype_exp$mic$tmm_voom[, ! colnames(celltype_exp$mic$tmm_voom) %in% sn_columbia$projid]
rownames(expr_mtx_sel) <- gsub("(*.)\\.(.*)","\\1", rownames(expr_mtx_sel))
expr_mtx_sel$ensembl = rownames(expr_mtx_sel)
gene_mod_map = data.frame(ensembl = rownames(expr_mtx_sel))
gene_mod_map = gene_mod_map %>% inner_join(unique(net_output[,c("ensembl","module_clusters")]), by = c("ensembl"))
gene_mod_map$cluster_lv3 = gsub("(.*)_M(.*)","\\2",gene_mod_map$module_clusters)

expr_mtx_sel_ord = expr_mtx_sel[match(gene_mod_map$ensembl,rownames(expr_mtx_sel)),] # order the matrix
identical(rownames(expr_mtx_sel_ord), gene_mod_map$ensembl) # must be TRUE 
## [1] TRUE
colors_mod = gene_mod_map$cluster_lv3
expr_mtx_sel_ord$ensembl = NULL
expr_mtx_sel_ord_t <- data.frame(t(as.matrix(expr_mtx_sel_ord)))
external_ME = moduleEigengenes(expr_mtx_sel_ord_t, colors = colors_mod)
mod_average = external_ME$averageExpr # data.frame

# Save results 
save(external_ME, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleEigengenes_lv3_mic.Rdata"))
write.table(external_ME$eigengenes, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleEigengenes_lv3_mic.txt"), sep = "\t", quote = F, row.names = T)
write.table(external_ME$averageExpr, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleaverageExpr_lv3_mic.txt"), sep = "\t", quote = F, row.names = T)

Regressions

data4linear_reg <- mod_average # match codes name 
phenotype_dt = phenotypes[match(rownames(data4linear_reg), phenotypes$projid), ]

all(rownames(data4linear_reg) == phenotype_dt$projid) # Must be TRUE. Match IDs
## [1] TRUE
# hist(phenotype_dt$amyloid_sqrt)
message(paste0("Unique participants: ", length(unique(phenotype_dt$projid)))) 
## Unique participants: 190
res_test = run_module_trait_association(data4linear_reg, phenotype_dt, pheno_list, covariates = c("age_death","msex", "educ")) # covariates to adjust
matrix_rsquared = res_test$matrix_rsquared
matrix_pvalue = res_test$matrix_pvalue
save(res_test, file = paste0(work_dir, "replication_analysis_filt/MIT_res_test_mic.Rdata"))

plot_module_trait_association_heatmap(res_test)

Top results

Top result by covariate.

# res_test$all_stats_df 
createDT(res_test$all_stats_df %>% group_by(phenotype) %>% slice_head(n = 1))

Nominal pvalue

createDT(matrix_pvalue)

Significant results

Threshold: At lemic one module with adjusted pvalue < 0.05.

plot_module_trait_association_heatmap(res_test, show_only_significant = T, signif_cutoff = 0.05)

Excitatory neurons

net_output = unique(sn_modules[sn_modules$cell_type == "ext", ])
expr_mtx_sel = celltype_exp$ext$tmm_voom[, ! colnames(celltype_exp$ext$tmm_voom) %in% sn_columbia$projid]
rownames(expr_mtx_sel) <- gsub("(*.)\\.(.*)","\\1", rownames(expr_mtx_sel))
expr_mtx_sel$ensembl = rownames(expr_mtx_sel)
gene_mod_map = data.frame(ensembl = rownames(expr_mtx_sel))
gene_mod_map = gene_mod_map %>% inner_join(unique(net_output[,c("ensembl","module_clusters")]), by = c("ensembl"))
gene_mod_map$cluster_lv3 = gsub("(.*)_M(.*)","\\2",gene_mod_map$module_clusters)

expr_mtx_sel_ord = expr_mtx_sel[match(gene_mod_map$ensembl,rownames(expr_mtx_sel)),] # order the matrix
identical(rownames(expr_mtx_sel_ord), gene_mod_map$ensembl) # must be TRUE 
## [1] TRUE
colors_mod = gene_mod_map$cluster_lv3
expr_mtx_sel_ord$ensembl = NULL
expr_mtx_sel_ord_t <- data.frame(t(as.matrix(expr_mtx_sel_ord)))
external_ME = moduleEigengenes(expr_mtx_sel_ord_t, colors = colors_mod)
mod_average = external_ME$averageExpr # data.frame

# Save results 
save(external_ME, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleEigengenes_lv3_ext.Rdata"))
write.table(external_ME$eigengenes, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleEigengenes_lv3_ext.txt"), sep = "\t", quote = F, row.names = T)
write.table(external_ME$averageExpr, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleaverageExpr_lv3_ext.txt"), sep = "\t", quote = F, row.names = T)

Regressions

data4linear_reg <- mod_average # match codes name 
phenotype_dt = phenotypes[match(rownames(data4linear_reg), phenotypes$projid), ]

all(rownames(data4linear_reg) == phenotype_dt$projid) # Must be TRUE. Match IDs
## [1] TRUE
# hist(phenotype_dt$amyloid_sqrt)
message(paste0("Unique participants: ", length(unique(phenotype_dt$projid)))) 
## Unique participants: 191
res_test = run_module_trait_association(data4linear_reg, phenotype_dt, pheno_list, covariates = c("age_death","msex", "educ")) # covariates to adjust
matrix_rsquared = res_test$matrix_rsquared
matrix_pvalue = res_test$matrix_pvalue
save(res_test, file = paste0(work_dir, "replication_analysis_filt/MIT_res_test_ext.Rdata"))

plot_module_trait_association_heatmap(res_test)

Top results

Top result by covariate.

# res_test$all_stats_df 
createDT(res_test$all_stats_df %>% group_by(phenotype) %>% slice_head(n = 1))

Nominal pvalue

createDT(matrix_pvalue)

Significant results

Threshold: At least one module with adjusted pvalue < 0.05.

plot_module_trait_association_heatmap(res_test, show_only_significant = T, signif_cutoff = 0.05)

Inhibitory neurons

net_output = unique(sn_modules[sn_modules$cell_type == "inh", ])
expr_mtx_sel = celltype_exp$inh$tmm_voom[, ! colnames(celltype_exp$inh$tmm_voom) %in% sn_columbia$projid]
rownames(expr_mtx_sel) <- gsub("(*.)\\.(.*)","\\1", rownames(expr_mtx_sel))
expr_mtx_sel$ensembl = rownames(expr_mtx_sel)
gene_mod_map = data.frame(ensembl = rownames(expr_mtx_sel))
gene_mod_map = gene_mod_map %>% inner_join(unique(net_output[,c("ensembl","module_clusters")]), by = c("ensembl"))
gene_mod_map$cluster_lv3 = gsub("(.*)_M(.*)","\\2",gene_mod_map$module_clusters)

expr_mtx_sel_ord = expr_mtx_sel[match(gene_mod_map$ensembl,rownames(expr_mtx_sel)),] # order the matrix
identical(rownames(expr_mtx_sel_ord), gene_mod_map$ensembl) # must be TRUE 
## [1] TRUE
colors_mod = gene_mod_map$cluster_lv3
expr_mtx_sel_ord$ensembl = NULL
expr_mtx_sel_ord_t <- data.frame(t(as.matrix(expr_mtx_sel_ord)))
external_ME = moduleEigengenes(expr_mtx_sel_ord_t, colors = colors_mod)
mod_average = external_ME$averageExpr # data.frame

# Save results 
save(external_ME, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleEigengenes_lv3_inh.Rdata"))
write.table(external_ME$eigengenes, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleEigengenes_lv3_inh.txt"), sep = "\t", quote = F, row.names = T)
write.table(external_ME$averageExpr, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleaverageExpr_lv3_inh.txt"), sep = "\t", quote = F, row.names = T)

Regressions

data4linear_reg <- mod_average # match codes name 
phenotype_dt = phenotypes[match(rownames(data4linear_reg), phenotypes$projid), ]

all(rownames(data4linear_reg) == phenotype_dt$projid) # Must be TRUE. Match IDs
## [1] TRUE
# hist(phenotype_dt$amyloid_sqrt)
message(paste0("Unique participants: ", length(unique(phenotype_dt$projid)))) 
## Unique participants: 190
res_test = run_module_trait_association(data4linear_reg, phenotype_dt, pheno_list, covariates = c("age_death","msex", "educ")) # covariates to adjust
matrix_rsquared = res_test$matrix_rsquared
matrix_pvalue = res_test$matrix_pvalue
save(res_test, file = paste0(work_dir, "replication_analysis_filt/MIT_res_test_inh.Rdata"))

plot_module_trait_association_heatmap(res_test)

Top results

Top result by covariate.

# res_test$all_stats_df 
createDT(res_test$all_stats_df %>% group_by(phenotype) %>% slice_head(n = 1))

Nominal pvalue

createDT(matrix_pvalue)

Significant results

Threshold: At least one module with adjusted pvalue < 0.05.

plot_module_trait_association_heatmap(res_test, show_only_significant = T, signif_cutoff = 0.05)

Oligodendrocytes

net_output = unique(sn_modules[sn_modules$cell_type == "oli", ])
expr_mtx_sel = celltype_exp$oli$tmm_voom[, ! colnames(celltype_exp$oli$tmm_voom) %in% sn_columbia$projid]
rownames(expr_mtx_sel) <- gsub("(*.)\\.(.*)","\\1", rownames(expr_mtx_sel))
expr_mtx_sel$ensembl = rownames(expr_mtx_sel)
gene_mod_map = data.frame(ensembl = rownames(expr_mtx_sel))
gene_mod_map = gene_mod_map %>% inner_join(unique(net_output[,c("ensembl","module_clusters")]), by = c("ensembl"))
gene_mod_map$cluster_lv3 = gsub("(.*)_M(.*)","\\2",gene_mod_map$module_clusters)

expr_mtx_sel_ord = expr_mtx_sel[match(gene_mod_map$ensembl,rownames(expr_mtx_sel)),] # order the matrix
identical(rownames(expr_mtx_sel_ord), gene_mod_map$ensembl) # must be TRUE 
## [1] TRUE
colors_mod = gene_mod_map$cluster_lv3
expr_mtx_sel_ord$ensembl = NULL
expr_mtx_sel_ord_t <- data.frame(t(as.matrix(expr_mtx_sel_ord)))
external_ME = moduleEigengenes(expr_mtx_sel_ord_t, colors = colors_mod)
mod_average = external_ME$averageExpr # data.frame

# Save results 
save(external_ME, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleEigengenes_lv3_oli.Rdata"))
write.table(external_ME$eigengenes, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleEigengenes_lv3_oli.txt"), sep = "\t", quote = F, row.names = T)
write.table(external_ME$averageExpr, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleaverageExpr_lv3_oli.txt"), sep = "\t", quote = F, row.names = T)

Regressions

data4linear_reg <- mod_average # match codes name 
phenotype_dt = phenotypes[match(rownames(data4linear_reg), phenotypes$projid), ]

all(rownames(data4linear_reg) == phenotype_dt$projid) # Must be TRUE. Match IDs
## [1] TRUE
# hist(phenotype_dt$amyloid_sqrt)
message(paste0("Unique participants: ", length(unique(phenotype_dt$projid)))) 
## Unique participants: 191
res_test = run_module_trait_association(data4linear_reg, phenotype_dt, pheno_list, covariates = c("age_death","msex", "educ")) # covariates to adjust
matrix_rsquared = res_test$matrix_rsquared
matrix_pvalue = res_test$matrix_pvalue
save(res_test, file = paste0(work_dir, "replication_analysis_filt/MIT_res_test_oli.Rdata"))

plot_module_trait_association_heatmap(res_test)

Top results

Top result by covariate.

# res_test$all_stats_df 
createDT(res_test$all_stats_df %>% group_by(phenotype) %>% slice_head(n = 1))

Nominal pvalue

createDT(matrix_pvalue)

Significant results

Threshold: At least one module with adjusted pvalue < 0.05.

plot_module_trait_association_heatmap(res_test, show_only_significant = T, signif_cutoff = 0.05)

Endothelial cells

net_output = unique(sn_modules[sn_modules$cell_type == "end", ])
expr_mtx_sel = celltype_exp$end$tmm_voom[, ! colnames(celltype_exp$end$tmm_voom) %in% sn_columbia$projid]
rownames(expr_mtx_sel) <- gsub("(*.)\\.(.*)","\\1", rownames(expr_mtx_sel))
expr_mtx_sel$ensembl = rownames(expr_mtx_sel)
gene_mod_map = data.frame(ensembl = rownames(expr_mtx_sel))
gene_mod_map = gene_mod_map %>% inner_join(unique(net_output[,c("ensembl","module_clusters")]), by = c("ensembl"))
gene_mod_map$cluster_lv3 = gsub("(.*)_M(.*)","\\2",gene_mod_map$module_clusters)

expr_mtx_sel_ord = expr_mtx_sel[match(gene_mod_map$ensembl,rownames(expr_mtx_sel)),] # order the matrix
identical(rownames(expr_mtx_sel_ord), gene_mod_map$ensembl) # must be TRUE 
## [1] TRUE
colors_mod = gene_mod_map$cluster_lv3
expr_mtx_sel_ord$ensembl = NULL
expr_mtx_sel_ord_t <- data.frame(t(as.matrix(expr_mtx_sel_ord)))
external_ME = moduleEigengenes(expr_mtx_sel_ord_t, colors = colors_mod)
mod_average = external_ME$averageExpr # data.frame

# Save results 
save(external_ME, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleEigengenes_lv3_end.Rdata"))
write.table(external_ME$eigengenes, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleEigengenes_lv3_end.txt"), sep = "\t", quote = F, row.names = T)
write.table(external_ME$averageExpr, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleaverageExpr_lv3_end.txt"), sep = "\t", quote = F, row.names = T)

Regressions

data4linear_reg <- mod_average # match codes name 
phenotype_dt = phenotypes[match(rownames(data4linear_reg), phenotypes$projid), ]

all(rownames(data4linear_reg) == phenotype_dt$projid) # Must be TRUE. Match IDs
## [1] TRUE
# hist(phenotype_dt$amyloid_sqrt)
message(paste0("Unique participants: ", length(unique(phenotype_dt$projid)))) 
## Unique participants: 179
res_test = run_module_trait_association(data4linear_reg, phenotype_dt, pheno_list, covariates = c("age_death","msex", "educ")) # covariates to adjust
matrix_rsquared = res_test$matrix_rsquared
matrix_pvalue = res_test$matrix_pvalue
save(res_test, file = paste0(work_dir, "replication_analysis_filt/MIT_res_test_end.Rdata"))

plot_module_trait_association_heatmap(res_test)

Top results

Top result by covariate.

# res_test$all_stats_df 
createDT(res_test$all_stats_df %>% group_by(phenotype) %>% slice_head(n = 1))

Nominal pvalue

createDT(matrix_pvalue)

Significant results

Threshold: At least one module with adjusted pvalue < 0.05.

plot_module_trait_association_heatmap(res_test, show_only_significant = T, signif_cutoff = 0.05)

OPCs

net_output = unique(sn_modules[sn_modules$cell_type == "opc", ])
expr_mtx_sel = celltype_exp$opc$tmm_voom[, ! colnames(celltype_exp$opc$tmm_voom) %in% sn_columbia$projid]
rownames(expr_mtx_sel) <- gsub("(*.)\\.(.*)","\\1", rownames(expr_mtx_sel))
expr_mtx_sel$ensembl = rownames(expr_mtx_sel)
gene_mod_map = data.frame(ensembl = rownames(expr_mtx_sel))
gene_mod_map = gene_mod_map %>% inner_join(unique(net_output[,c("ensembl","module_clusters")]), by = c("ensembl"))
gene_mod_map$cluster_lv3 = gsub("(.*)_M(.*)","\\2",gene_mod_map$module_clusters)

expr_mtx_sel_ord = expr_mtx_sel[match(gene_mod_map$ensembl,rownames(expr_mtx_sel)),] # order the matrix
identical(rownames(expr_mtx_sel_ord), gene_mod_map$ensembl) # must be TRUE 
## [1] TRUE
colors_mod = gene_mod_map$cluster_lv3
expr_mtx_sel_ord$ensembl = NULL
expr_mtx_sel_ord_t <- data.frame(t(as.matrix(expr_mtx_sel_ord)))
external_ME = moduleEigengenes(expr_mtx_sel_ord_t, colors = colors_mod)
mod_average = external_ME$averageExpr # data.frame

# Save results 
save(external_ME, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleEigengenes_lv3_opc.Rdata"))
write.table(external_ME$eigengenes, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleEigengenes_lv3_opc.txt"), sep = "\t", quote = F, row.names = T)
write.table(external_ME$averageExpr, file = paste0(work_dir, "replication_analysis_filt/MIT_moduleaverageExpr_lv3_opc.txt"), sep = "\t", quote = F, row.names = T)

Regressions

data4linear_reg <- mod_average # match codes name 
phenotype_dt = phenotypes[match(rownames(data4linear_reg), phenotypes$projid), ]

all(rownames(data4linear_reg) == phenotype_dt$projid) # Must be TRUE. Match IDs
## [1] TRUE
# hist(phenotype_dt$amyloid_sqrt)
message(paste0("Unique participants: ", length(unique(phenotype_dt$projid)))) 
## Unique participants: 191
res_test = run_module_trait_association(data4linear_reg, phenotype_dt, pheno_list, covariates = c("age_death","msex", "educ")) # covariates to adjust
matrix_rsquared = res_test$matrix_rsquared
matrix_pvalue = res_test$matrix_pvalue
save(res_test, file = paste0(work_dir, "replication_analysis_filt/MIT_res_test_opc.Rdata"))

plot_module_trait_association_heatmap(res_test)

Top results

Top result by covariate.

# res_test$all_stats_df 
createDT(res_test$all_stats_df %>% group_by(phenotype) %>% slice_head(n = 1))

Nominal pvalue

createDT(matrix_pvalue)

Significant results

Threshold: At least one module with adjusted pvalue < 0.05.

plot_module_trait_association_heatmap(res_test, show_only_significant = T, signif_cutoff = 0.05)

Session info

sessionInfo()
## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Stream 8
## 
## Matrix products: default
## BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] RColorBrewer_1.1-3    circlize_0.4.16       ComplexHeatmap_2.15.4
##  [4] performance_0.10.0    lmerTest_3.1-3        lme4_1.1-28          
##  [7] Matrix_1.6-1.1        WGCNA_1.70-3          fastcluster_1.2.3    
## [10] dynamicTreeCut_1.63-1 janitor_2.1.0         furrr_0.3.1          
## [13] future_1.33.2         kableExtra_1.3.4      lubridate_1.9.3      
## [16] forcats_1.0.0         stringr_1.5.1         dplyr_1.1.4          
## [19] purrr_1.0.2           readr_2.1.5           tidyr_1.3.1          
## [22] tibble_3.2.1          ggplot2_3.5.0         tidyverse_2.0.0      
## [25] Seurat_5.0.0          SeuratObject_5.0.0    sp_2.1-1             
## 
## loaded via a namespace (and not attached):
##   [1] utf8_1.2.4             spatstat.explore_3.2-5 reticulate_1.35.0     
##   [4] tidyselect_1.2.1       RSQLite_2.3.3          AnnotationDbi_1.56.2  
##   [7] htmlwidgets_1.6.4      Rtsne_0.16             munsell_0.5.1         
##  [10] preprocessCore_1.56.0  codetools_0.2-18       ica_1.0-3             
##  [13] DT_0.30                miniUI_0.1.1.1         withr_3.0.0           
##  [16] spatstat.random_3.2-1  colorspace_2.1-0       progressr_0.14.0      
##  [19] Biobase_2.54.0         highr_0.10             knitr_1.46            
##  [22] rstudioapi_0.16.0      stats4_4.1.2           ROCR_1.0-11           
##  [25] tensor_1.5             listenv_0.9.1          GenomeInfoDbData_1.2.7
##  [28] polyclip_1.10-6        bit64_4.0.5            parallelly_1.37.1     
##  [31] vctrs_0.6.5            generics_0.1.3         xfun_0.43             
##  [34] timechange_0.3.0       R6_2.5.1               doParallel_1.0.17     
##  [37] GenomeInfoDb_1.30.1    clue_0.3-65            bitops_1.0-7          
##  [40] spatstat.utils_3.0-4   cachem_1.0.8           promises_1.3.0        
##  [43] scales_1.3.0           nnet_7.3-16            gtable_0.3.4          
##  [46] Cairo_1.6-0            globals_0.16.3         processx_3.8.4        
##  [49] goftest_1.2-3          spam_2.10-0            rlang_1.1.3           
##  [52] systemfonts_1.0.6      GlobalOptions_0.1.2    splines_4.1.2         
##  [55] lazyeval_0.2.2         impute_1.68.0          checkmate_2.0.0       
##  [58] spatstat.geom_3.2-7    yaml_2.3.8             reshape2_1.4.4        
##  [61] abind_1.4-5            crosstalk_1.2.1        backports_1.4.1       
##  [64] httpuv_1.6.15          Hmisc_4.6-0            tools_4.1.2           
##  [67] jquerylib_0.1.4        BiocGenerics_0.40.0    ggridges_0.5.4        
##  [70] Rcpp_1.0.12            plyr_1.8.9             base64enc_0.1-3       
##  [73] zlibbioc_1.40.0        RCurl_1.98-1.13        ps_1.7.6              
##  [76] rpart_4.1-15           deldir_1.0-9           GetoptLong_1.0.5      
##  [79] pbapply_1.7-2          cowplot_1.1.3          S4Vectors_0.32.4      
##  [82] zoo_1.8-12             ggrepel_0.9.5          cluster_2.1.2         
##  [85] magrittr_2.0.3         data.table_1.15.4      RSpectra_0.16-1       
##  [88] scattermore_1.2        lmtest_0.9-40          RANN_2.6.1            
##  [91] fitdistrplus_1.1-11    matrixStats_1.3.0      hms_1.1.3             
##  [94] patchwork_1.2.0        mime_0.12              evaluate_0.23         
##  [97] xtable_1.8-4           jpeg_0.1-9             shape_1.4.6.1         
## [100] fastDummies_1.7.3      IRanges_2.28.0         gridExtra_2.3         
## [103] compiler_4.1.2         KernSmooth_2.23-20     crayon_1.5.2          
## [106] websocket_1.4.1        minqa_1.2.6            htmltools_0.5.8.1     
## [109] later_1.3.2            tzdb_0.4.0             Formula_1.2-4         
## [112] DBI_1.2.2              MASS_7.3-54            boot_1.3-28           
## [115] cli_3.6.2              insight_0.18.5         parallel_4.1.2        
## [118] dotCall64_1.1-0        igraph_2.0.3           pkgconfig_2.0.3       
## [121] numDeriv_2016.8-1.1    foreign_0.8-81         plotly_4.10.4         
## [124] spatstat.sparse_3.0-3  xml2_1.3.6             foreach_1.5.2         
## [127] svglite_2.1.3          bslib_0.7.0            chromote_0.2.0        
## [130] XVector_0.34.0         webshot_0.5.2          rvest_1.0.4           
## [133] snakecase_0.11.0       digest_0.6.35          sctransform_0.4.1     
## [136] RcppAnnoy_0.0.21       spatstat.data_3.0-3    Biostrings_2.62.0     
## [139] rmarkdown_2.26         leiden_0.4.3           htmlTable_2.4.0       
## [142] uwot_0.1.16            shiny_1.8.1.1          rjson_0.2.21          
## [145] nloptr_2.0.3           lifecycle_1.0.4        nlme_3.1-153          
## [148] jsonlite_1.8.8         viridisLite_0.4.2      fansi_1.0.6           
## [151] pillar_1.9.0           lattice_0.20-45        KEGGREST_1.34.0       
## [154] fastmap_1.1.1          httr_1.4.7             survival_3.2-13       
## [157] GO.db_3.14.0           glue_1.7.0             png_0.1-8             
## [160] iterators_1.0.14       bit_4.0.5              stringi_1.8.3         
## [163] sass_0.4.9             blob_1.2.4             RcppHNSW_0.5.0        
## [166] latticeExtra_0.6-29    memoise_2.0.1          irlba_2.3.5.1         
## [169] future.apply_1.11.2